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Abstract 

When a strong magnetostatic field is present, vacuum effectively appears as a linear, uni- 
axial, dielectric-magnetic medium for small-magnitude optical fields. The availability of the 
frequency-domain dyadic Green function when the magnetostatic field is spatially uniform 
facilitates the formulation of an integral equation for the scattering of an optical field by a 
spatially varying magnetostatic field in vacuum. This integral equation can be numerically 
treated by using the method of moments as well as the coupled dipole method. Furthermore, 
the principle underlying the strong-property-fluctuation theory allows the homogenization of 
a spatially varying magnetostatic field in the context of light scattering. 

Key words: coupled dipole method, depolarization, homogenization, method of moments, 
quantum electrodynamics. 
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1 Introduction 



In classical electrodynamics, light propagating in vacuum (i.e., matter free space) is not considered 
to be affected by the presence of a magnetostatic field. This is because classical vacuum is a linear 
medium wherein the principle of superposition holds. But, in quantum electrodynamics (QED), 
vacuum is a nonlinear medium (Jackson, 1998). It can, however, be linearized for a rapidly time- 
varying electromagnetic field with a small amplitude in the presence of a slowly varying (or static) 
magnetic field (Adler, 2007). The price of linearization is that the QED vacuum appears as an 
anisotropic dielectric-magnetic medium for optical fields (Adler, 1971). 

Our modest aim in this paper is to derive an integral equation for the scattering of a high- 
frequency electromagnetic field (typified by light) by a strong magnetostatic field in vacuum. For 
this purpose, we exploit the analytical machinery developed during the last two decades for the 
frequency-domain analysis of electromagnetic fields in complex mediums (Singh & Lakhtakia, 
2000; Weiglhofer & Lakhtakia, 2003). As the derived integral equation shall have to be solved 
numerically in general, we also outline the solution strategies provided by the method of moments 
(Miller, Medgyesi-Mitschang, & Newman, 1991; Wang, 1991) and the coupled dipole method 
(Purcell & Pennypacker, 1973; Lakhtakia, 1990). Finally, we adopt the principle underlying the 
strong property-fluctuation theory to homogenize the magnetostatic field in the context of light 
scattering. 

A note about notation: 3-vectors (6-vectors) are in normal (bold) face and underlined, whereas 
3x3 clyadics (6x6 dyadics) are in normal (bold) face and double underlined. The position vector 
is denoted by r = xu^ + yUy + zu^ in a Cartesian coordinate system with unit vectors u^, Uy, and 
u^, whereas time is denoted by t. The real part of a complex-valued quantity ( is written as 3? {(}. 

2 Constitutive Equations for Optical Fields 

Suppose that all space is matter-free and that a magnetostatic field B^'^{r) is present everywhere. 
Optical fields (superscripted "o" in this paper) satisfy the source-free Maxwell equations (Adler, 
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1971) 



V-S^fci) = 
V-^°(r,i) = 

Vx^°(r,i) = -^S°(r,t) 



(1) 



VxH°{r,t) = -§-^D'>{r,t) 

provided \B°\ «; |B'''=(r)| Vr. As per the linearized version of QED, the optical fields appearing in 
the foregoing equations obey the constitutive equations (Adler,1971) 



(2) 



D''{r,t)=eog''ir)-E''{r,t) 

H"{r,t) = fx^'g''{r)-B"{r,t) 

where eo = 8.8542 x 10~^^ F and /xq = 47r x 10"'^ H are constants characterizing the 
classical vacuum, whereas the relative permittivity dyadic 



e°(r) = 1 - 8eo4^B'^''{r) ■ B^^{r) / + 28eQ6f,^ B^\r)B^\r) 



dc( 



-ydc( 



and the relative impermeablity dyadic^ 



:°(r) = 1 - 8eo4^B^\r) ■ B'^^r) I_ - l(Seo4^ B^''{r)B^''{r) 



(3) 



(4) 



emerge from QED to dictate the influence of the magnetostatic field of high magnitude on the 
optical field. Here, cq = 1/ ^/eo^o is the speed of light in classical vacuum, whereas 



- = 8.3229 X 10-32 i^g-i g2 



(5) 



^ 45(47reo)2c7 

contains the electronic charge e^e = 1.6022 x 10"^^ C, the electronic mass niei = 9.1096 x 10"^^ kg, 
and the reduced Planck constant h = 1.0546 x 10-^* J s. Clearly, the QED vacuum appears to the 
optical field as a spatiotemporally local, spatially nonhomogoncous, temporally unvarying, uniaxial 
dielectric- magnetic medium. Set ft = to convert from the QED vacuum to the classical vacuum, 
and the influence of the magnetostatic field on the optical field vanishes. 



Impermeability is the reciprocal of the permeability. 
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3 Dyadic Green Function for Uniform Magnetostatic Field 



Suppose that Bf'^{r) = is spatially uniform everywhere, so that the QED vacuum then is 
spatially homogeneous. Defining 

ac = eocU\Bf 



we can write the relative permittivity dyadic 



(6) 



1^ = (1 - Sac){I - u^u^) + (1 + 20ac) 



and the relative permeability dyadic 



24a. 



(7) 



(8) 



for the QED vacuum, with / as the 3x3 identity dyadic. The propagation of optical plane waves, 
characterized by 

E"{r, i) = 3? exp [i{k-r- Lot)] } , (9) 

etc., with amplitude vector £^°, wave vector fc, and angular frequency w was analyzed by Adler 
(1971) using standard mathematical techniques. Both the relative permittivity and the relative 
permeability dyadics are uniaxial and share the same distinguished axis. Therefore, the medium is 
optically birefringent except when k is parallel to u.; furthermore, both plane waves propagating 
in a fixed direction are to be classified as extraordinary in optical parlance (Lakhtakia, Varadan, 
& Varadan, 1991). 

For a more general consideration of frequency-domain optical fields, it is best to write 



E°{r,t) = 5R{e°(r)exp(-ia;f)} , 



(10) 



etc., where e°(r) is a phasor. The Maxwell curl equations may then be compactly recast with the 
help of 6-vectors and 6x6 dyadics as 



L(V)+ia;C° • T (r) = s°(r) , 



where 



L(V) = 



2 V X I 

-V X / 



(11) 



(12) 
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c° 



(13) 



(14) 



s°(r) represents the (high-frequency) electric and magnetic source current densities that engender 
f°(r), and is the 3x3 null dyadic. 

Equation (11) is converted into the integral equation 



(15) 



where f/^(r) is the solution of (11) when the source term on its right side is null-valued everywhere. 
The 6x6 dyadic Green function g°(r — r') is the solution of 



L(V) + iujC" ■ g°{r, r') =IS{r- r') , 



(16) 



where I is the 6x6 identity dyadic and S{r — r'} is Dirac delta function. 

The 6x6 dyadic Green fimction g"(r, r') is available from the literature on electromagnetic fields 
in uniaxial mediums as (Weiglhofer, 1990) 

L -(iweoip-^ • (V x I) 

(zw/uo/x")"^ • (V X X) I 
where the 3x3 dyadic Green functions 



f{r,r') = 



9" ir,r') 

=ce ~ 

g° ir,r') 



, (17) 



9" {r,r') = -{l-8ac)-^{{iweo)-' fvV + /c^(l + 20ac)(a-'l Tefcr') + ^^Mo ¥'(r:,r')| (18) 

=ce L L — J = J 



and 



9° (r,r') = -(l-8ac){(ia;Mo)-' [vV + fc^(l-24ae)-^(/x°)-'l 7^(r,r') - iweo <^(r:,r')| (19) 

=cm L L =c J = J 

employ the usual wavenumber ko — w(eoMo)^^^ for classical vacuum. The 3x3 dyadic function 



1 + 20ac 
1 - 8oc 



7°(z:,z:')-(y 



1 - 8ac 



24o. 



7^(e,z:') 



{R X u^){R X 
l-R X u^P 



+ 



I - u^Uc - 2 



{R X Uc){R X 



and the scalar functions 



l-R X u^p 



iko\R X u^p 



exp(iA;o^^e) 
4nR° 



(20) 



(21) 
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and 



contain 



^ exp(zfcoj?^) 



R = r-r' , 



T30 _ 



1 + 20ac 
1 - 8ac 

1 - 8ac 
1 - 24a, 



1/2 



1/2 



(22) 

(23) 
(24) 

(25) 



4 Formulation of Integral Equation 

Returning to the more general case of a nonuniform magnetostatic field Bf'^{r), and making use 
of the representation (10), we see that the two curl equations in (1) may be written together as 



^(V)+ia;£°(r)] ■r{r) = 0, 



(26) 



where is the null 6-vector and 



(27) 



eo|°(r:) Q 

g MO [g"{r)r' 

As a simple Green-function formalism is unlikely to be found for this more general case, we define 
a suitable uniform magnetostatic field and proceed as follows. 



4.1 Derivation 

Equation (26) is recast as 

where 



L(V)+za;g° •r(r)=s°<,(r). 



(28) 



(29) 



is an equivalent high-frequency source current density 6-vector that contains the spatial variations 
of Bf'^(r) in relation to B^. Clearly, should be chosen carefully. If the variations of Bf'^{r) 
are confined to some bounded region only, should be chosen as the value of Bf'^{r) outside 
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that region. If B^^{r) varies over all space, may emerge from some homogenization or spatial- 
averaging procedure. Once that choice has been made, the solution of (28) — and therefore of 
(26) — may be written as 



i''{L)=il{r)+ioJ JdY£{r,r') 



C°-C°(r') •f°(r'), 



(30) 



where f^(ri) is the homogeneous part of the solution of (28). 

The integrand on the right side of (30) is singular at r = r' and therefore requires additional 
treatment. An ellipsoidal region containing the location r midway between its two focuses is 
identified in the domain of integration. The surface of Vg is the set of points 



r,{eg,(j)g) =r + eg-Uq, 61^ e [0,7r] , 0g e [0,27r] 
where e is a positive scalar, the unit vector 

Ug = (Uj. cos (pq + Uy sin (j)q) sin Og + cos 6q , 



(31) 



(32) 



and the dyadic 

^ = ttx UxlLx + y^yiky + W^M^ (33) 

has positive eigenvalues a^, ay, and az, all less than or equal to unity. The integral is divided into 
two parts: the first is to be evaluated over 14 only, the second over all space except V^. Both parts 
are evaluated in the limit e ^ in order to get the desired integral equation 

= fhiD + ^wD • \C° - C°(r)l • T (r) + PV / d\' g°(r, r') • fc° - C°(r')l • r{r') , (34) 



wherein the symbol PV identifies that the integral following it has to be evaluated in the 'principal- 
value sense', as described. 



4.2 Depolcirization Dyadic 

The 6x6 depolarization dyadic D in (34) is a two-dimensional surface integral (Michel & Weiglhofer 
1997; Weiglhofer, Lakhtakia, & Michel 1997): 

r2iT 



D = 



d4lq 



dOg sin^q 



(u-u-^){&-u-^) 



= (u-u-^yno-iu-u-^) ^ 



(35) 
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In general, it has to be evaluated numerically, which can accomplished quite easily by using the 
Gauss-Legendre quadrature scheme, for instance (Mackay & Weiglhofer, 2000). 

In two special situations, D can be obtained analytically. Suppose, first, that the ellipsoidal 
region is spherical. Then, C/ = 7 so that 



D = - 



27r 



MOST 











(36) 



For e° defined in (7) and /i° defined in (8), the two integrals on the right side of (36) can be 
evaluated in closed form to yield (Mackay & Lakhtakia 2005) 



where 



and 



D = - 



lUJ 

A-Kkl 







l-8ac 

eo(l - 8ac) 



I I I psinh" 



-1 /i-p- 
p 



2\l-p 

1 / psec^^ 1 



rc(p) = <^ 



2 V {p-iy p~^, 
1 



p 



(l_p)2 1-p 

1 sec~^ yfp 



for < p < 1 

for p > 1 
for < p < 1 

for p > 1 



_ l+20ae 



l-8ao 



1^ 1-24(1,. 



> . 



(37) 



(38) 



(39) 



(40) 



(41) 



The second special situation comes when one of the semi-axes of the ellipsoidal region 
coincides with u^. Without loss of generality, let us choose u^. = and then introduce the 3x3 
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dyadics 



!£"' = ur^ ■ !±° ■ ur^ = i4 ^Lxikx + iJ-i uyUy + 14 uzUz 



wherein 



l-8ac 



l-8ac 
5 — 



y a: 



,o' _ l+20ac 



^a; (l-8a<:)a2 



3/ ~ (l-8ac)aj 



> . 



I^z — (l-24oc)a^ - 

It then transpires that the depolarization dyadic may be represented as (Weiglhofer, 1998) 



eour^-dt-ur^ 



where 



with 



dZ = 



d7. 



Vy -Vx I Vz - Vy 



V°z -V"x 



1/2 



, i^(Ai, A2) - , E{X^,X2) 



Vy 



1/2 



S(Ai,A2) 



(42) 



(43) 



(44) 



(45) 



(46) 



which involve F(Ai, A2) and ^^(Ai, A2) as elliptic integrals of the first and second kinds (Gradshteyn 
& Ryzhik, 1980), respectively, with arguments 



Ai = tan ^ 



A2 



Vz - Vx 



m [Vy - Tlx 

%'{v"z-v°x': 



1/2 



1/2 



(47) 
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4.3 Numerical— Solution Techniques 

As may be easily guessed, (34) shall have to be solved numerically. When the spatial variations of 

Bf'^{r) are confined to some bounded region Vi whereas the magnctostatic field is uniform with a 
value everywhere outside Vi, the method of moments (Miller, Medgyesi Mitschang, & Newman, 
1991; Wang, 1991) and the coupled dipolc method (PurccU & Pcnnypackcr, 1973; Lakhtakia,1990) 
offer relatively easy algorithms to implement. Both methods are related to each other (Lakhtakia, 
1992), and their adaptations for the present purposes are described as follows. 



4.3.1 Method of Moments 

Implementation of the method of moments requires that Vi be replaced by a lattice of points r„, 
n G Attached to every r„ is an electrically small region described by a shape dyadic and 

of volume w„; the sum {J2n=i equals the volume of Vi. Furthermore, a depolarization dyadic 
is associated with every r„. 
Equation (34) is specialized to r = r„, n G [1, A''], to obtain the set of algebraic equations 

N 

+ E ^'-Cfen.rJ- [g-mzijj-rdj, ne[i,Ar], (48) 

where f^(r) represents the incident optical field (i.e., the optical fields when B^'^{r) = Vr e Vi). 
This set of equations is rewritten as 



N 



(49) 



m=l 



where 



qMOm ^ 1 J _ . \ci - C°(r jl } 5„„ - t;„g°(r„, r„) • [c° - C^rj] (1 - 5nra) , (50) 

L~ — L — ^ — J J =c L — ^ — J 

with 6nm as the Kronecker delta function. 

Equations (49) are solved for all f°(r„), n G [l,N]hy a variety of numerical techniques (Carna- 
han, Luther, & Wilkes, 1969), with the conjugate gradient method being the preferred technique 
when N is large (Strang, 1986; Sarkar, 1991). Once all £"(?:„) are known, the total field at any 
r^Vi may be determined from (34) as 

N 

r(r) ~ mr) + E ^rn£{r,rj ■ [q^ - Q^rj] • r{rj , r^Vi. (51) 

m=l 
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The difference f°(r) — f^(r) yields the scattered optical field at r ^ V^. 



4.3.2 Coupled Dipole Method 

In the coupled dipole method, the first and the third terms on the right side of (48) are used to 
define an exciting optical field via 

N 

fLc(z:n)^f^(rJ+ E ^m£{r„,rJ.[Ql-Q''{rj].r{rJ, nG[l,N]. (52) 
Thereafter, the set of equations (48) is rewritten as 



r(r:J«fLc(in) + ^'^D 



C°(rJ -rdlj, ne[l,iV], 



whence 



r (r J ~ {l - • [g - £°(r: j] } ' • fLcfcn) , n e [1, TV] . 



(53) 



(54) 



Clearly, ^^^Jyr_^ contains not only the incident optical field at r„ but also the scattered optical 
field due to all other lattice points in V^;. 
Substituting (54) in (52), we obtain 

AT 



f;:(rJ=E Q''''''-f-cfcJ, ne[l,iV] 



m=l 



where 



contains 



—=nm ~ — 



(55) 



(56) 



(57) 



Equations (55) have to be solved numerically too — just like (49) — in order to determine f excfcn)' 
n e [1, N]. Thereafter, by virtue of (51) and (54), the the total optical field at any r^Vi may be 
determined from (34) as 

N 



(58) 



The reason for the name of this technique is that the product v^X (56) and (58) may be consid- 
ered as a 6x6 polarizability dyadic of the electrically small region associated with (Weiglhofer, 
Lakhtakia, & Michel, 1997; Lakhtakia & Weiglhofer, 2000). 
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4.4 Homogenization of Magnetostatic Field 



The coupled dipole method naturally leads to the homogenization problem of replacing a spatially 
varying Bf'^{r) in some region V hy a uniform magnetostatic field 5^^^ for the analysis of optical 
fields in V. 

Following the strong property fluctuation theory (SPFT) (Tsang & Kong, 1981; Michel & 
Lakhtakia, 1995; Mackay, Lakhtakia, & Weiglhofer 2000), let us begin by introducing the con- 
tinuous analog of the exciting optical field of (54) as 



f Led) = {l - ■ - £°(r)] } • r (r) , reV, 



(59) 



where C° . is defined in analogy with C° — as per (7), (8), and (13) — and D is to be calculated 
by using C° = C°^. and [/ = / on the right side of (35). Thereby, the integral equation (34) may 
be rewritten as 

CciD = lUD + io^ PV /_ g" (r, r') • x{l) ■ fLcfcO , (60) 

JV =urn = 

where g° {LiI!) is defined in analogy with g°(r, r') in (16) with C° = C° ., and the polarizability 

=uni =c — ^ — um 

density dyadic is defined as 



C° .-C°(r)j • (l-iwD- [C° .-C°(r)]) 

^=.um = J 1^ = = \=.um = J J 



(61) 



in analogy to x ■ 

Within the SPFT framework, the volume integrals of both sides of (60) are used to compute the 
constitutive parameters of a uniform medium which approximate to those of the medium described 
by C°(r). The simplest estimate of these constitutive parameters is provided by C°^., which is 
found by imposing the condition 

!d\x{r)=g. (62) 

JV 

I : order to deduce B_^^^ from (62), we introduce the piecewise-uniform approximation 

Bf\r) = B^, r&Vj, (63) 
with V = Uj=i %■ The corresponding uniform 6x6 constitutive dyadics are defined as 



revs-, 



(64) 
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with 3x3 dyadic components 

1° = (1 - Seoc^e B^. . B^) I + 28eocl^B.Bj 

g = [(1 - SeocUR^ ■ B^) L - IGeocUB^B^] 
Then, the condition (62) impUes that 



(66) 

where fj is the volumetric proportion of region Vj relative to V. This nonlinear dyadic equation 
can be straightforwardly solved for . — and thereby — by standard numerical procedures 
(Michel, Lakhtakia, & Weiglhofcr, 1998). 

Higher-order corrections to C°^^ can also be derived using the SPFT, depending upon the 
statistical details of the spatial fluctuations in Bf'^{r) (Mackay, Lakhtakia, & Weiglhofer, 2000). 

5 Concluding Remarks 

To conchide, we have formulated an integral equation for the scattering of light (and other high- 
frequency) electromagnetic radiation in vacuum by either a strong magnetostatic field or a strong 
low-frequency magnetic field. This integral equation, set in the language of 6-vectors and 6x6 
dyadics, exploits the frequency-domain dyadic Green function for a linear, homogeneous, uniaxial, 
dielectric-magnetic medium. The derived integral equation can be solved numerically, using the 
method of moments as well as the coupled dipole method. Finally, we have shown that the strong- 
property-fluctuation theory allows the homogenization of a spatially varying magnetostatic field 
in the present context. 

Acknowledgement: AL thanks John Collins of the Physics Department, Penn State, for a dis- 
cussion. TGM is supported by a Royal Society of Edinburgh/ Scottish Executive Support Research 
Fellowship. 
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